Langerhans AUTO

#Starting from saved object

library(Seurat)
library(ggplot2)
library(Matrix)
library(dplyr)


markCluster <- function(fex){
  FeaturePlot(langer, features = fex, cols = c("grey", "red"), reduction = "tsne", combine = F, label = T)
  VlnPlot(langer, features = fex, split.by = "seurat_clusters", pt.size = 0, combine = F)
}

langer <- readRDS('/home/vqf/proyectos/TFGs/Ruben/NoteBooks/langer.rds')
DimPlot(langer, reduction = 'tsne', split.by =  'genotype', label = T, repel = T)

Get characteristic genes for each cluster. The tutorial advises using the original data.

TSNEPlot(langer)

Lmna expression

FeaturePlot(langer, features = c("Lmna"), cols = c("grey", "red"), reduction = "tsne", combine = F, label = T)

[[1]]

VlnPlot(langer, features = c("Lmna"), split.by = "seurat_clusters", pt.size = 0)

VlnPlot(langer, features = c("Lmna"), split.by = "genotype", pt.size = 0)

Calr expression

FeaturePlot(langer, features = c("Calr"), cols = c("grey", "red"), reduction = "tsne", combine = F, label = T)

[[1]]

VlnPlot(langer, features = c("Calr"), split.by = "seurat_clusters", pt.size = 0)

VlnPlot(langer, features = c("Calr"), split.by = "genotype", pt.size = 0)

Hspa5 expression

FeaturePlot(langer, features = c("Hspa5"), cols = c("grey", "red"), reduction = "tsne", combine = F, label = T)

[[1]]

VlnPlot(langer, features = c("Hspa5"), split.by = "seurat_clusters", pt.size = 0)

VlnPlot(langer, features = c("Hspa5"), split.by = "genotype", pt.size = 0)

Ins expression

fex <- c("Ins1", "Ins2", "Gcg", "Ppy", "Sst", "Rbp4", "Mafb", "Mafa", "Ghrl" , "Zfp608", "Irx2", "Nr4a2")
DotPlot(langer, features = c("Gcg", "Ins1", "Sst", "Ppy"), group.by = "seurat_clusters")

markCluster(fex)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

Proliferating

fex <- c("Gmnn", "Mki67", "Pcna")
markCluster(fex)

[[1]]

[[2]]

[[3]]

Dying

fex <- c("Casp3", "Bcl2", "Parp1", "Nfkb1")
markCluster(fex)

[[1]]

[[2]]

[[3]]

[[4]]

Exocrine

fex <- c("Amy2a1", "Cela3a", "Cpa1")
markCluster(fex)

[[1]]

Endocrine

fex <- c("Chga", "Chgb", "Pcsk1","Pcsk2", "Scg2", "Scg5", "Neurod1", "Isl1", "Pax6")
markCluster(fex)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

Fibroblast

fex <- c("Col3a1", "Fbn1", "Lum")
markCluster(fex)

[[1]]

[[2]]

[[3]]

Smooth-muscle cells

fex <- c("Acta2", "Actg2", "Myh11")
markCluster(fex)

[[1]]

[[2]]

[[3]]

SMC and Fibroblast tend to be in mixed clusters so we can select more acurrate threshold

RidgePlot(subset(langer, seurat_clusters==9), features = c("Acta2"))+ geom_vline(xintercept=0.5)

RidgePlot(subset(langer, seurat_clusters==9), features = c("Lum"))+ geom_vline(xintercept=0.37)

FeatureScatter(subset(langer, seurat_clusters==9), "Acta2", "Lum", span = T)+ geom_vline(xintercept=0.5)+ geom_hline(yintercept=0.37)

RidgePlot(subset(langer, seurat_clusters==9), features = c("Acta2"))+ geom_vline(xintercept=0.7)

RidgePlot(subset(langer, seurat_clusters==9), features = c("Col3a1"))+ geom_vline(xintercept=0.45)

FeatureScatter(subset(langer, seurat_clusters==9), "Acta2", "Col3a1", span = T)+ geom_vline(xintercept=0.5)+ geom_hline(yintercept=0.7)

Endothelial cells

fex <- c("Cd34", "Pecam1", "Sele", "Slc2a1", "Vwf")
markCluster(fex)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

Leukocytes

fex <- c("Ptprc")
markCluster(fex)

[[1]]

B cells

fex <- c("Cd19", "Cr2", "Ms4a1")
markCluster(fex)

[[1]]

[[2]]

[[3]]

T cells

fex <- c("Cd3e", "Cd4", "Cd8a", "Foxp3", "L17a", "Gzmb", "Tox", "Pdcd1", "Slamf6")
markCluster(fex)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

Plasma cells

fex <- c("Pim2")
markCluster(fex)

[[1]]

Dendritic cells

fex <- c("Gzmb", "Il3ra")
markCluster(fex)

[[1]]

[[2]]

Schwann cells

fex <- c("Mbp", "Mpz", "S100b", "Sox10")
markCluster(fex)

[[1]]

[[2]]

Monocytes

fex <- c("Fcgr4")
markCluster(fex)

[[1]]

NK cells

fex <- c("Cd27")
markCluster(fex)

[[1]]

Macrophages

fex <- c("Ccr5", "Cd68", "Marco", "Mrc1", "Msr1")
markCluster(fex)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

Ductal

fex <- c("Krt17")
markCluster(fex)

[[1]]

Beta and delta cells

fex <- c("Ins1", "Sst")
markCluster(fex)

[[1]]

[[2]]

We select ours clusters with Beta and Delta cells

RidgePlot(subset(langer, seurat_clusters==20), features = c("Ins1"))+ geom_vline(xintercept=5.5)

RidgePlot(subset(langer, seurat_clusters==20), features = c("Sst"))+ geom_vline(xintercept=4)

FeatureScatter(subset(langer, seurat_clusters==20), "Ins1", "Sst", span = T)+ geom_vline(xintercept=5.5)+ geom_hline(yintercept=4)

Alpha and PP

fex <- c("Gcg", "Ppy")
markCluster(fex)

[[1]]

[[2]]

We select our cluster with Alpha and PP cells

RidgePlot(subset(langer, seurat_clusters==23), features = c("Gcg"))+ geom_vline(xintercept=5.5)

RidgePlot(subset(langer, seurat_clusters==23), features = c("Ppy"))+ geom_vline(xintercept=2.5)

FeatureScatter(subset(langer, seurat_clusters==23), "Gcg", "Ppy", span = T)+ geom_vline(xintercept=5.5)+ geom_hline(yintercept=2.5)

FeatureScatter(subset(langer, seurat_clusters==23), "Gcg", "Sst", span = T)+ geom_vline(xintercept=5.5)+ geom_hline(yintercept=2.7)

##Clasification

Identify cell types

  d <- GetAssayData(langer)
  m <- matrix(data=c("Other"),nrow=ncol(d),ncol=1,byrow = F,dimnames = list(colnames(d),c("celltype")))
  setClusterName <- function(seuratObj, mobj, clusters, clusterName){
    ind <- seuratObj@meta.data[seuratObj@meta.data$seurat_clusters %in% clusters,]
    mobj[rownames(ind), "celltype"] <- c(clusterName)
    return(mobj)
  }

Umbrales

# NK cells
ind <- subset(langer, subset = Cd27 > 0.2 )
m[rownames(ind@meta.data), "celltype"] <- "NK cell"

# Plasma cells
ind <- subset(langer, subset = Pim2 > 0.2 )
m[rownames(ind@meta.data), "celltype"] <- "Plasma cell"

# Macrophages
ind <- subset(langer, subset = Cd68 > 0.3  | Marco > 0.3| Mrc1 > 0.3 | Msr1 > 0.3 )
m[rownames(ind@meta.data), "celltype"] <- "Macrophage"

# Smooth muscle cells
ind <- subset(langer, subset = (Acta2 > 0.5 & Lum > 0.37) | (Acta2 > 0.7 & Col3a1 > 0.45))
m[rownames(ind@meta.data), "celltype"] <- "SMC"

# Fibroblasts
ind <- subset(langer, subset = (Acta2 < 0.5 & Lum < 0.37) | (Acta2 < 0.7 & Col3a1 < 0.45))
m[rownames(ind@meta.data), "celltype"] <- "Fibroblast"

# Leukocytes
ind <- subset(langer, subset = Ptprc > 1)
m[rownames(ind@meta.data), "celltype"] <- "Leukocyte"

# T cells
ind <- subset(langer, subset = Cd3e > 0.3 | Foxp3 > 0.3 | Cd8a > 0.3 | Cd4 > 0.3 )
m[rownames(ind@meta.data), "celltype"] <- "T cell"

# B cells
ind <- subset(langer, subset = Ms4a1 > 0.3 | Cd19 > 0.3)
m[rownames(ind@meta.data), "celltype"] <- "B cell"

# Endothelial cells
ind <- subset(langer, subset = Cd34 > 0.3 | Slc2a1 > 0.3 | Pecam1 > 0.3)
m[rownames(ind@meta.data), "celltype"] <- "Endothelial"

# Beta cells
ind <- subset(langer, subset = (Ins1 > 5.5 & Gcg < 2 & Ppy < 2))
m[rownames(ind@meta.data), "celltype"] <- "Beta"

# Alpha cells
ind <- subset(langer, subset = (Gcg > 5.5 & Ppy < 2.5 & Ins1 < 5.5))
m[rownames(ind@meta.data), "celltype"] <- "Alpha"

# Delta cells
ind <- subset(langer, subset = (Sst > 4 & Ins1 < 5.5 & Gcg < 2))
m[rownames(ind@meta.data), "celltype"] <- "Delta"

# PP cells
ind <- subset(langer, subset = (Ppy > 2.5 & Gcg < 5.5 & Ins1 < 2))
m[rownames(ind@meta.data), "celltype"] <- "PP"

# BD cells
ind <- subset(langer, subset = (Ins1 > 5.5 & Sst > 4 & Ppy < 2))
m[rownames(ind@meta.data), "celltype"] <- "BD"
  
# AP cells
ind <- subset(langer, subset = (Gcg > 4 & Ppy > 2 & Ins1 < 1.5 & Sst < 1.5))
m[rownames(ind@meta.data), "celltype"] <- "AP"
  

langer <- AddMetaData(langer, m, col.name = "celltype")
Idents(langer) <- 'celltype'
TSNEPlot(langer, label=T, split.by = 'genotype')

Save islets cells.

langer$celltype <- m[rownames(langer@meta.data), "celltype"]
islets <- subset(langer, celltype=='Alpha' | celltype=='Beta' | celltype=='Delta' | celltype=="DB" | celltype=="BD" | celltype== "AP"| celltype== "PP" | celltype== "PA")
saveRDS(islets, file = '/home/vqf/proyectos/TFGs/Ruben/NoteBooks/islets.rds')

Interpretate Data

table(langer@meta.data$genotype)

   KO    WT 
31588 26893 

tbl_wt <- table(langer@meta.data$celltype[langer@meta.data$genotype == "WT"])
cbind(tbl_wt,round(prop.table(tbl_wt)*100,2))

            tbl_wt      
Alpha           64  0.24
AP              26  0.10
B cell        3692 13.73
BD               2  0.01
Beta            68  0.25
Delta           15  0.06
Endothelial   4373 16.26
Fibroblast    5455 20.28
Leukocyte     1377  5.12
Macrophage       2  0.01
NK cell          1  0.00
Other          224  0.83
PP              45  0.17
SMC            101  0.38
T cell       11448 42.57

tbl_ko <- table(langer@meta.data$celltype[langer@meta.data$genotype == "KO"])
cbind(tbl_ko,round(prop.table(tbl_ko)*100,2))

            tbl_ko      
Alpha          215  0.68
AP             144  0.46
B cell        3384 10.71
BD               2  0.01
Beta           266  0.84
Delta           36  0.11
Endothelial   5084 16.09
Fibroblast    6037 19.11
Leukocyte     1356  4.29
Macrophage       4  0.01
Other          409  1.29
Plasma cell      2  0.01
PP             164  0.52
SMC            293  0.93
T cell       14192 44.93

tbl <- table(islets@meta.data$celltype, islets@meta.data$orig.ident)
cbind(tbl,round(prop.table(tbl)*100,2))

      158WT 161KO 235KO 238WT 243KO 244WT 158WT 161KO 235KO 238WT 243KO 244WT
Alpha    33    64    94    26    57     5  3.15  6.11  8.98  2.48  5.44  0.48
AP       10    35    82    12    27     4  0.96  3.34  7.83  1.15  2.58  0.38
BD        1     0     1     1     1     0  0.10  0.00  0.10  0.10  0.10  0.00
Beta     18    53   113    22   100    28  1.72  5.06 10.79  2.10  9.55  2.67
Delta    10    16     9     5    11     0  0.96  1.53  0.86  0.48  1.05  0.00
PP       27    59    67    15    38     3  2.58  5.64  6.40  1.43  3.63  0.29

RidgePlot(subset(langer, seurat_clusters==20), features = c("Ins1"))+ geom_vline(xintercept=5.5)

RidgePlot(subset(langer, seurat_clusters==20), features = c("Sst"))+ geom_vline(xintercept=4)

Inmuno <- subset(langer, celltype== 'Macrophage' | celltype== 'Leukocyte' | celltype== 'B cell' | celltype== 'T cell'| celltype== 'NK cell')
DimPlot(islets, reduction = "tsne", split.by = "genotype", group.by = "celltype", label = TRUE, pt.size = 1.4) + ggtitle("Distribución de células inmunes (tSNE)")

DimPlot(langer, reduction = "tsne", group.by = "celltype", label = TRUE, pt.size = 1.2) +
  ggtitle("Distribución de tipos celulares (tSNE)")

# DimPlot sólo con células endocrinas
p <-DimPlot(islets, reduction = "tsne", split.by = "genotype", group.by = "celltype", label = TRUE, pt.size = 1.4) +
  ggtitle("Distribución de células endocrinas (tSNE)")

p + xlim (20, 30) + ylim(15, 35)